/*
 * @(#)BroydenFletcherGoldfarbShanno.java        1.0 2000/06/26
 *
 * This file is part of Xfuzzy 3.0, a design environment for fuzzy logic
 * based systems.
 *
 * (c) 2000 IMSE-CNM. The authors may be contacted by the email address:
 *                    xfuzzy-team@imse.cnm.es
 *
 * Xfuzzy is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by
 * the Free Software Foundation.
 *
 * Xfuzzy is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 * for more details.
 */

package xfuzzy.xfsl.model.algorithm;

import xfuzzy.xfsl.model.*;
import xfuzzy.lang.*;
import xfuzzy.xfds.XfdsDataSet;

/**
 * Algoritmo de Broyden, Fletcher, Goldfarb y Shanno
 * 
 * @author Francisco Jos� Moreno Velo
 * 
 * @see "Fletcher, R., Practical Methods of Optimization, 
 * John Wiley & Sons, Ltd., 1986."
 *
 */
public class BroydenFletcherGoldfarbShanno extends XfslGradientBasedAlgorithm {

	//----------------------------------------------------------------------------//
	//                            MIEMBROS PRIVADOS                               //
	//----------------------------------------------------------------------------//

	/**
	 * Estimaci�n del hessiano
	 */
	private double[][] lh;
	
	/**
	 * Tolerancia en la b�squeda lineal
	 */
	private double tol;
	
	/**
	 * L�mite de iteraciones para la b�squeda lineal
	 */
	private int limit;

	//----------------------------------------------------------------------------//
	//                                CONSTRUCTOR                                 //
	//----------------------------------------------------------------------------//

	/**
	 * Constructor
	 */
	public BroydenFletcherGoldfarbShanno() {
		this.tol = 0.1;
		this.limit = 10;
	}

	//----------------------------------------------------------------------------//
	//                             M�TODOS P�BLICOS                               //
	//----------------------------------------------------------------------------//

	/**
	 * Devuelve el c�digo de identificaci�n del algoritmo
	 */
	public int getCode() {
		return BFGS;
	}

	/**
	 * Actualiza los par�metros de configuraci�n del algoritmo
	 */
	public void setParameters(double[] param) throws XflException {
		if(param.length != 2) throw new XflException(26);
		tol = test(param[0], REDUCE);
		limit = (int) test(param[1], INTEGER);
	}

	/**
	 * Obtiene los par�metros de configuraci�n del algoritmo 
	 */
	public XfslAlgorithmParam[] getParams() {
		XfslAlgorithmParam[] pp = new XfslAlgorithmParam[2];
		pp[0] = new XfslAlgorithmParam(tol, 0.1, REDUCE, "Line-search Tolerance");
		pp[1] = new XfslAlgorithmParam(limit, 10, INTEGER, "Search Iteration Limit");
		return pp;
	}
	
	/**
	 * Ejecuta una iteraci�n del algoritmo
	 */
	public XfslEvaluation iteration(Specification spec, XfdsDataSet pattern,
			XfslErrorFunction ef) throws XflException {
		XfslEvaluation prev = computeErrorGradient(spec,pattern,ef);
		OptimizingFunction function = new OptimizingFunction(spec,pattern,ef);
		Parameter[] param = spec.getAdjustable();
		double[] pt = new double[param.length];
		for(int i=0; i<param.length; i++) pt[i] = param[i].value;
		double[][] h;

		if(init) { init=false; h=newHessian(param); }
		else h = BFGS(param);

		double[] g = new double[param.length];
		for(int i=0; i<param.length; i++) g[i] = -param[i].getDeriv();
		double[] p = product(h,g);

		XfslEvaluation eval = function.linmin(p,prev,tol,limit);

		for(int i=0; i<param.length; i++) {
			param[i].setPrevDesp(param[i].value - pt[i]);
			param[i].setPrevDeriv(-g[i]);
			param[i].setDeriv(0);
		}

		lh = h;
		return eval;
	}

	//----------------------------------------------------------------------------//
	//                             M�TODOS PRIVADOS                               //
	//----------------------------------------------------------------------------//

	/**
	 * Inicializa el valor del Hessiano
	 */
	private double[][] newHessian(Parameter[] param) {
		double[][] h = new double[param.length][param.length];
		for(int i=0; i<param.length; i++) h[i][i] = 1.0;
		return h;
	}

	/**
	 * Multiplica el Hessiano por el gradiente	
	 */
	private double[] product(double[][] h, double[] g) {
		double[] p = new double[h.length];
		for(int i=0; i<h.length; i++)
			for(int j=0; j<g.length; j++)
				p[i] += h[i][j]*g[j];
		return p;
	}

	/**
	 * Actualizacion de Broyden, Fletcher, Goldfarb y Shanno
	 */
	 private double[][] BFGS(Parameter[] param) {
		 double[] dx = new double[param.length];
		 double[] dg = new double[param.length];
		 double xg=0;
		 for(int i=0; i<param.length; i++) {
			 dx[i] = param[i].getPrevDesp();
			 dg[i] = param[i].getDeriv() - param[i].getPrevDeriv();
			 xg += dx[i]*dg[i];
		 }
		 if(xg == 0) return newHessian(param);

		 double[] hg = product(lh,dg);
		 double ghg=0;
		 for(int i=0; i<param.length; i++) ghg += dg[i]*hg[i];

		 double alpha = 1+ghg/xg;
		 double[][] h = new double[param.length][param.length];
		 for(int i=0; i<param.length; i++)
			 for(int j=0; j<param.length; j++)
				 h[i][j] = lh[i][j] + (alpha*dx[i]*dx[j] - dx[i]*hg[j] - dx[j]*hg[i])/xg;
		 return h;
	 }

}

